library(pacman)
pacman::p_load(tidyverse, data.table, broom,hrbrthemes, plyr,latex2exp,openxlsx,ggpubr,ggpmisc)
rm(list=ls())
###############################################

###Spatial units
df_su = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "spatial_units"))[, list(spatial_id_name)]

###Prices
df = read.csv(file = paste0('results_ic/lf/p_i_c.csv'), header=FALSE)
setnames(df, old=c('V1','V2'), new=c('beef','crops') )
df_lf_P = df

df = read.csv(file = paste0('results_ic/dt/p_i_c.csv'), header=FALSE)
setnames(df, old=c('V1','V2'), new=c('beef','crops') )
df_cf_P = df

df = 100*(df_cf_P-df_lf_P)/df_lf_P
df$origin=df_su$spatial_id
dP_pp = gather(df, commodity, delta, beef:crops, factor_key=TRUE)

df = df_lf_P
df$origin=df_su$spatial_id
lf_P = gather(df, commodity, delta, beef:crops, factor_key=TRUE)

###Quantities
df = read.csv(file = paste0('results_ic/lf/Q.csv'), header=FALSE)
setnames(df, old=c('V1','V2'), new=c('beef','crops') )
df_lf_Q = df

df = read.csv(file = paste0('results_ic/dt/Q.csv'), header=FALSE)
setnames(df, old=c('V1','V2'), new=c('beef','crops') )
df_cf_Q = df

###Revenue
df_lf_R = df_lf_P*df_lf_Q
df_cf_R = df_cf_P*df_cf_Q
df = 100*(df_cf_R-df_lf_R)/df_lf_R
df$origin=df_su$spatial_id
dR_pp = gather(df, commodity, delta, beef:crops, factor_key=TRUE)

n_r_count = 1
for(n_r in c('farm-gate prices','farmer income')){
  
  if(n_r=='farm-gate prices'){df_rate = dP_pp}
  if(n_r=='farmer income'){df_rate = dR_pp}
  df_level = lf_P
  
  df = merge(df_level,df_rate,by=c('origin','commodity'))
  setnames(df, old=c('delta.x','delta.y'), new=c('level','rate'))
  df$income_definition= n_r 
  df=setDT(df)
  
  df$percentile = 100*rank(df$level)/length(df$level)
  df[commodity=='beef',]$percentile= 100*rank(df[commodity=='beef',]$level)/length(df[commodity=='beef',]$level)
  df[commodity=='crops',]$percentile= 100*rank(df[commodity=='crops',]$level)/length(df[commodity=='crops',]$level)

  if(n_r_count==1){df_r=df}else{df_r=rbind(df,df_r)}
  n_r_count=n_r_count+1
}

###Scatter plots
df = df_r[commodity=='beef',]
df$y_var=df$rate
df$x_var=df$percentile
df_full=df

df=df_full[income_definition %in% c('farm-gate prices'),]
ylab=TeX('$\\Delta$ farm-gate prices (pp)')
xlab = 'baseline income (percentile)'
min_val = quantile(df$x_var,0.025)
max_val = quantile(df$x_var,0.975)
df=df[x_var>min_val & x_var<max_val,]
fig_farmgate_prices=ggplot(data = df, aes(x = x_var, y=y_var) ) +
  stat_smooth(method = "lm", formula = y ~ x + I(x^2) , linewidth = 1, se=T, color='black')+
  stat_summary_bin(fun='mean', color='black', geom='point', alpha=0.5, size=2)+
  theme_minimal() +
  labs(x=xlab, y=ylab)+
  theme(axis.title.x = element_text(size=8, hjust = 1), 
        axis.title.y = element_text(size=8, hjust = 0.95),
        plot.title = element_text(size = 8.5, hjust = 0.5, vjust=2),
        axis.text.x = element_text(size=5),
        axis.text.y = element_text(size=5),
        aspect.ratio = 1,
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.position = 'bottom',
        legend.text=element_text(size=9))+
  guides(size = "none")

df=df_full[income_definition %in% c('farmer income'),]
ylab=TeX('$\\Delta$ farmer income (pp)')
xlab = 'baseline income (percentile)'
min_val = quantile(df$x_var,0.025)
max_val = quantile(df$x_var,0.975)
df=df[x_var>min_val & x_var<max_val,]
fig_farmer_income=ggplot(data = df, aes(x = x_var, y=y_var) ) +
  stat_smooth(method = "lm", formula = y ~ x + I(x^2) , linewidth = 1, se=T, color='black')+
  stat_summary_bin(fun='mean', color='black', geom='point', alpha=0.5, size=2)+
  theme_minimal() +
  labs(x=xlab, y=ylab)+
  theme(axis.title.x = element_text(size=8, hjust = 1), 
        axis.title.y = element_text(size=8, hjust = 0.95),
        plot.title = element_text(size = 8.5, hjust = 0.5, vjust=2),
        axis.text.x = element_text(size=5),
        axis.text.y = element_text(size=5),
        aspect.ratio = 1,
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.position = 'bottom',
        legend.text=element_text(size=9))+
  guides(size = "none")

figure=ggarrange(fig_farmgate_prices,fig_farmer_income, nrow=1, ncol=2)
ggsave(paste0("../../output/figures/appendix/figure6_regressivity.pdf"), width = 5, height = 2.5)
